Readme

The following script provides the general framework used to analyse demographic buffering continuum in Santos et al. (in review). Population responses to environmental stochasticity are primarily driven by survival-reproduction trade-offs and mediated by aridity.

Intermediary code and data steps are sourced along the framework.
Intermediary code and data include:
- 1 - Data cleaning and selection and its intermediary data 1 - Data cleaning and selection.R
- 2 - Life history traits calculation and its intermediary data 2 - Life history traits calculation.R
- 3 - Climatic variables extraction and estimation of stochastic elasticity
- ChelsaData Download and extraction - Google Colab.ipynb
- 3 - Climatic variables calculation.R
- 4 - Core function: Stochastic elasticities of variance
- MainFunction - Stochastic elasticities of variance lower level.R
- 4 - Model selection and MCMCglmm inference
- 5a - MCMCglmm - Model selection.R and 5a - MCMCglmm - Model selection.ipynb
- 5b- MCMCglmm run.R and 5b- MCMCglmm run.ipynb

1 - Setting

0.1. SessionInfo and package versions

Show all relevant information for R and package versions

## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=Portuguese_Brazil.utf8  LC_CTYPE=Portuguese_Brazil.utf8   
## [3] LC_MONETARY=Portuguese_Brazil.utf8 LC_NUMERIC=C                      
## [5] LC_TIME=Portuguese_Brazil.utf8    
## 
## time zone: America/Sao_Paulo
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] MCMCglmm_2.36     ape_5.8           coda_0.19-4.1     Matrix_1.7-0     
##  [5] ggtree_3.14.0     treeio_1.30.0     tidybayes_3.0.7   rstatix_0.7.2    
##  [9] vegan_2.6-8       lattice_0.22-6    permute_0.9-7     Rage_1.6.0       
## [13] Rcompadre_1.3.0   factoextra_1.0.7  FactoMineR_2.11   viridis_0.6.5    
## [17] viridisLite_0.4.2 scales_1.3.0      popdemo_1.3-1     popbio_2.8       
## [21] lubridate_1.9.3   forcats_1.0.0     stringr_1.5.1     dplyr_1.1.4      
## [25] purrr_1.0.2       readr_2.1.5       tidyr_1.3.1       tibble_3.2.1     
## [29] ggplot2_3.5.1     tidyverse_2.0.0  
## 
## loaded via a namespace (and not attached):
##  [1] gridExtra_2.3        sandwich_3.1-1       rlang_1.1.4         
##  [4] magrittr_2.0.3       multcomp_1.4-26      compiler_4.4.1      
##  [7] mgcv_1.9-1           vctrs_0.6.5          pkgconfig_2.0.3     
## [10] arrayhelpers_1.1-0   fastmap_1.2.0        backports_1.5.0     
## [13] utf8_1.2.4           rmarkdown_2.29       tzdb_0.4.0          
## [16] xfun_0.49            cachem_1.1.0         aplot_0.2.3         
## [19] jsonlite_1.8.8       flashClust_1.01-2    broom_1.0.7         
## [22] parallel_4.4.1       cluster_2.1.6        R6_2.5.1            
## [25] bslib_0.8.0          stringi_1.8.4        RColorBrewer_1.1-3  
## [28] car_3.1-3            jquerylib_0.1.4      estimability_1.5.1  
## [31] Rcpp_1.0.13          knitr_1.49           zoo_1.8-12          
## [34] splines_4.4.1        timechange_0.3.0     tidyselect_1.2.1    
## [37] rstudioapi_0.17.1    abind_1.4-8          yaml_2.3.10         
## [40] codetools_0.2-20     withr_3.0.2          posterior_1.6.0     
## [43] evaluate_1.0.1       gridGraphics_0.5-1   survival_3.6-4      
## [46] ggdist_3.3.2         pillar_1.9.0         tensorA_0.36.2.1    
## [49] carData_3.0-5        DiagrammeR_1.0.11    checkmate_2.3.2     
## [52] DT_0.33              ggfun_0.1.7          distributional_0.5.0
## [55] generics_0.1.3       hms_1.1.3            tidytree_0.4.6      
## [58] munsell_0.5.1        xtable_1.8-4         leaps_3.2           
## [61] glue_1.7.0           cubature_2.1.1       lazyeval_0.2.2      
## [64] emmeans_1.10.5       scatterplot3d_0.3-44 tools_4.4.1         
## [67] fs_1.6.5             visNetwork_2.1.2     mvtnorm_1.2-5       
## [70] grid_4.4.1           colorspace_2.1-1     patchwork_1.2.0     
## [73] nlme_3.1-165         Formula_1.2-5        cli_3.6.3           
## [76] fansi_1.0.6          expm_0.999-9         svUnit_1.0.6        
## [79] corpcor_1.6.10       gtable_0.3.6         yulab.utils_0.1.8   
## [82] sass_0.4.9           digest_0.6.36        ggplotify_0.1.2     
## [85] ggrepel_0.9.5        TH.data_1.1-2        htmlwidgets_1.6.4   
## [88] htmltools_0.5.8.1    lifecycle_1.0.4      multcompView_0.1-10 
## [91] MASS_7.3-61

2. Data selection and cleaning - load population models & phylogenetic information

Script avaliable in “1 - Data cleaning and selection.R”

file.edit("1 - Data cleaning and selection.R")

Produces two datasets:
- CleanData.RDS:
  - Filter matrix singularity,
  - Select only matrices with fecundity,
  - Individual matrices only
- supertree.RDS:
  - Extract supertree from the MOSAIC database Bernard et al. (2023)

3. Life history traits (Calculate life history traits)

Script available in “2 - Life history traits calculation.R”

Produces the dataset LHtraits.RDS by:
- Calculate life history traits and detect outliers
  - Life history traits include:
    - D: Reproductive window (\(L\))
      - Temporal spread of reproduction throughout the life span as quantified by the Demetrius (1974) ’s entropy (\(S\)).
        
        High/low \(S\) values correspond to iteroparous/semelparous populations
    - La: Mean age at first reproduction (\(L_a\))
      - Average time from birth to reproduction
    - LaProp: Probability that an individual becomes reproductive along its lifetime (\(p_a\))
    - e: Mean life expectancy (\(\eta _e\))
    - growth: individual development (\(\gamma\))
    - H: Distribution of mortality risk along the life cycle (\(H\))
  - Outlier detection was done using Mahalanobis distance
    - mahal.dist The mahalanobis distance
    - is.outlier: TRUE or FALSE to indicate outliers detected with the Mahalanobis distance

3.1 Fast-slow continuum / Life history PCA

Perform PCA with life-history data;
Check eigenvalues;
Plot eigenvalues and correlogram.

##                                      comp 1    comp 2     comp 3     comp 4
## eigenvalue                         2.696769  1.173172  0.7162735  0.6095803
## percentage of variance            44.946155 19.552872 11.9378909 10.1596722
## cumulative percentage of variance 44.946155 64.499028 76.4369184 86.5965906
##                                       comp 5      comp 6
## eigenvalue                         0.4286934   0.3755112
## percentage of variance             7.1448894   6.2585200
## cumulative percentage of variance 93.7414800 100.0000000

## $Dim.1
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##        correlation      p.value
## H        0.7782491 6.313930e-59
## LaProb   0.7487164 2.802252e-52
## Growth   0.7104679 6.192170e-45
## D        0.6479642 3.320243e-35
## e        0.1633963 5.779568e-03
## La      -0.7610525 6.144154e-55
## 
## $Dim.2
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##        correlation       p.value
## e        0.9300967 1.118102e-124
## LaProb   0.3165349  4.990971e-08
## La       0.2377485  5.185214e-05
## H        0.1668020  4.826793e-03
## D       -0.2444222  3.123859e-05
## Growth  -0.2526037  1.645514e-05
## 
## $Dim.3
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##        correlation      p.value
## D        0.6282915 1.345055e-32
## La       0.3476003 1.735429e-09
## LaProb   0.2631009 6.998948e-06
## H       -0.3471587 1.825150e-09

## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

4. Climate variables and environmental PCA

Extract and store climatic data from CHELSAcruts (Karger et al. 2017) in an intermediary dataset, remove collinear variable, perform environmental PCA, and plot the PCA axes across the globe.

4.1. Extract and store climate information

The extracted climatic data is available in Data/Climate_df.RDS

Because data download and information extraction are quite time consuming, a dedicated Jupyter notebook to run on Google Colab is provided: ChelsaData Download and extraction - Google Colab.ipynb

Run this Google Colab by typing in your browser: https://githubtocolab.com/Ecosantos/Demogbuff-pops

4.1. Analyse climatic timeseries, calculate environmental stochasticity, and summarise information into a PCA

3 - Climatic variables calculation.R use the already extracted climatic data to:

Produces data summary and extract metrics described in methods
- “Mean_trend_TMax”¹:
- “Stoch_noisesize_TMax”
- “Ampli_trend_TMax”
- “Ampli_trend_TMin”
- “Mean_trend_Prec”
- “Ampli_trend_Prec”
Produces a new dedicated dataset: climate_df.RDS

4.1.1. Checking collinearity and retaining non-redundant climatic variables

4.1.2. Climatic/environmental PCA

Run the environmental PCA and show its output

##                                      comp 1    comp 2    comp 3     comp 4
## eigenvalue                         1.799101  1.231347  1.049691  0.8588607
## percentage of variance            30.077274 20.585597 17.548683 14.3583897
## cumulative percentage of variance 30.077274 50.662871 68.211554 82.5699438
##                                       comp 5     comp 6
## eigenvalue                         0.7230174   0.319578
## percentage of variance            12.0873684   5.342688
## cumulative percentage of variance 94.6573121 100.000000

## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

## $Dim.1
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##                      correlation       p.value
## Stoch_noisesize_TMax   0.8668087 6.183297e-100
## Ampli_trend_TMin       0.4590744  2.130222e-18
## Ampli_trend_Prec       0.2177142  7.386559e-05
## Ampli_trend_TMax       0.1981180  3.193055e-04
## Mean_trend_Prec       -0.4508056  1.008298e-17
## Mean_trend_TMax       -0.7434070  1.507656e-58
## 
## $Dim.2
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##                  correlation      p.value
## Ampli_trend_Prec   0.7886580 2.021982e-70
## Mean_trend_TMax    0.4022397 4.147529e-14
## Ampli_trend_TMin  -0.2481373 5.784008e-06
## Ampli_trend_TMax  -0.2935288 6.707701e-08
## Mean_trend_Prec   -0.5476851 6.527546e-27
## 
## $Dim.3
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##                  correlation      p.value
## Ampli_trend_TMax   0.7230280 5.494684e-54
## Ampli_trend_TMin   0.4402050 6.962149e-17
## Ampli_trend_Prec   0.3842174 6.557083e-13
## Mean_trend_TMax    0.3647761 1.071806e-11
## Mean_trend_Prec    0.2271603 3.473356e-05

4.3 Mapping the axes of the environmental PCA

Plot the spatial distribution of studies with available climatic data, as well as the values of the two main principal component analyses: life-history and environmental PCAs.

5. Calculate Demographic buffering from Vital Rates
(Lower Level Parameters)

5.1. Sourcing main function used in this calculation

The main function calculates the stochastic elasticities within respect to variance of lower-level vital rates

This script creates the following functions:
1. My.vitalRatePerturbation: calculates the stochastic elasticity with respect to variance \(| \sum E_v^\sigma |\)
  - This function is an adaptation of the VitalRatePerturbation function in the Rage package (see Jones et al. (2022))
2. array_to_matrix: An ancillary function to transform timeseries of MPMs, usually returned as an array element to a timeserie of MPM in an object class matrix.

## [1] "An ancillary function - 'array_to_matrix' - has been created, be pround!"
## [1] "Core function 'my.vitalRatePerturbation' - has been created, You nailed it!"
## 'my.vitalRatePerturbation' is based on: 
## - VitalRatePerturbation function in the Rage package: https://github.com/jonesor/Rage. 
## - Stochastic elasticticity in Caswell 2001 p.406-407. 
## Outputs are: 
## [[1]] Lower-level Stochastic Sensitivities 
## [[2]] Lower-level Stochastic Elasticities 
## [[3]] Lower-level Stochastic Elasticities 
## [[4]] Lower-level Stochastic Elasticities with respect  to variance
## [[5]] Matrix level stochastic Elasticitity with respect to the meanNULL

5.2. Derivative demographic buffering

Derivates demographic buffering in plants and animals
Merge into a single data frame

5.2.1. Number of populations with timeseries longer than three

## [1] 182

5.2.2. Merging Demographic Buffering into a single data frame

5.2.3. Produce a definitive demographic buffering data frame

This new data frame will be ready to merge with climatic data and life history traits

## Warning in cor(.): o desvio padrão é zero

6. Merging all (Climatic data + Life history + Buffering)

7. Building the super tree and making sure that works on phylogenetic analyses

Check must return:
- No duplicates: Each species must appear only once as a tip in the phylogeny. Duplicate taxa are not allowed.
- No polytomies: The phylogeny must be fully resolved, meaning no nodes with more than two descendants. Methods require bifurcating trees, so polytomies must be removed.
- is.rooted == TRUE: The tree must be rooted, indicating a clear ancestral lineage and direction of evolution. This defines an initial state in evolutionary time.
- is.binary == TRUE: All internal nodes must have exactly two descendants. This ensures the tree is strictly bifurcating, which is required for many comparative methods.
- is.ultrametric == TRUE: The tree must be ultrametric, meaning all tips are equidistant from the root. Branch lengths should reflect evolutionary time or history.
- any edge == 0 == FALSE: No branch (edge) in the tree should have a length of zero. Zero-length branches can cause computational problems or indicate unresolved relationships.

## [1] TRUE

## [1] FALSE

## ***************************************************************
## *                          Note:                              *
## *    force.ultrametric does not include a formal method to    *
## *    ultrametricize a tree & should only be used to coerce    *
## *   a phylogeny that fails is.ultrametric due to rounding --  *
## *    not as a substitute for formal rate-smoothing methods.   *
## ***************************************************************

## [1] FALSE

## [1] TRUE

## [1] TRUE

## [1] TRUE

## [1] TRUE

## [1] FALSE

## numeric(0)

8. Produce a final dataset and phylogenetic analyses

## [1] "Total pops"

## [1] 134

## [1] "Total species"

## [1] 89

## Kingdom
## Animalia  Plantae 
##       13      121

## Kingdom
## Animalia  Plantae 
##       11       78

## [1] "Total MPMs"

## [1] 889

Because animal population are scarce, we keep the analyses for plants only.

8.1. Prepare the associated phylogenetic information

9. GLMM parameters & data input

9.1. Prepare a separate dataset `data_model`

Data_model includes absolute values of sum of stochastic elasticity within respect to variance \(| \sum E_v^\sigma |\) (i.e, cumulativa_SigEla).²

save(data_model,subtree_Plants,
     file = "Data/GLMMdata.Rdata")

9.2. Selecting the best model

To optimise the model selection analysis a dedicated Jupyter notebook integrated with Google Colab was created: 5b - MCMCglmm run.ipynb

9.3. Determine the fixed effect of all models and response variables

Because we run multiple models representing the combination of several attributes (vital rates, ex. survival and reproduction) for the same fixed model, each element below is important to automatise the process. Once each element is created separately, we can run a generic model in a loop model(variable ~ fi

9.4. OPTIMISING GLMM MODELS

Because we run several models with and without phylogenetic correction, it takes a long time to run in a in-house computer using R with the script 5b - MCMCglmm run.R. Thus, to accelerate the model performance and save time, we created a jupyter notebook integrated with Google Colab 5b - MCMCglmm run.ipynb.

Most importantly, a partial output was created with the GLMM models in Data/MCMCglmm_output.rds.

The next steps use this partial output to continue running the analyses.

10. Buffering patterns

10.1 Ranking buffering

Rank more and less buffered populations. Again, the sum of stochastic elasticity within respect to variance \(| \sum E_v^\sigma |\) is used.

Plot

An accessory function is necessary to better represent very small numbers.

Plot

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

10.3. Proportional contribution

Quantify the proportional contribution of each vital rate to \(| \sum E_v^\sigma |\). We used percentage RelContrib*100 to improve readability in the main manuscript.

10.3.1 Proportional contribution by vital rate

10.3.1 Proportional contribution by life form

10.3.3. Plots

Raw plots are produced to be edited with proper software.

Plot Elasticities

The plot shows the stochastic elasticities of each vital rate (A) and its relative effect on \(\lambda_{s}\), it means \(E_{v}^\sigma\)

Plot proportional contribution

Proportional contribution of each \(|E_{v}^\sigma|\) to \(|\sum E_{v}^\sigma|\)

Raw figure

Figure 2 main manuscript

11. MCMCglmm results and tidy

Use Data/MCMCglmm_output.rds to check estimates and hypotheses testing

11.1 Loading GLMM outputs

## $Simple_models
## [1] "Plants"
## 
## $Phylogenetic_models
## [1] "Plants"

11.2. Data Harmonisation

11.2.1 As a summary

Use this for a data summary along the main manuscript

## <list_of<
##   tbl_df<
##     Trait     : character
##     Taxa      : character
##     Model     : character
##     Statistics: character
##     post.mean : double
##     low95     : double
##     high95    : double
##     eff.samp  : double
##     pMCMC     : double
##     sig       : character
##   >
## >[5]>
## [[1]]
## # A tibble: 12 × 10
##    Trait  Taxa  Model Statistics post.mean    low95  high95 eff.samp pMCMC sig  
##    <chr>  <chr> <chr> <chr>          <dbl>    <dbl>   <dbl>    <dbl> <dbl> <chr>
##  1 Clona… Plan… Phylo LHPC.1     -0.00263  -0.00871 0.00338   10581. 0.393 Non-…
##  2 Clona… Plan… Phylo LHPC.2     -0.00101  -0.0110  0.00946    9605. 0.853 Non-…
##  3 Clona… Plan… Phylo ClimPC.1    0.00152  -0.0103  0.0127     9900. 0.794 Non-…
##  4 Clona… Plan… Phylo ClimPC.2   -0.00396  -0.0153  0.00831    9900. 0.509 Non-…
##  5 Clona… Plan… Phylo LHPC.1:LH…  0.00123  -0.00421 0.00662    9900  0.655 Non-…
##  6 Clona… Plan… Phylo ClimPC.1:… -0.00170  -0.00930 0.00572    8527. 0.656 Non-…
##  7 Clona… Plan… Simp… LHPC.1     -0.000547 -0.00740 0.00607    9900  0.875 Non-…
##  8 Clona… Plan… Simp… LHPC.2     -0.00591  -0.0175  0.00682    9900. 0.337 Non-…
##  9 Clona… Plan… Simp… ClimPC.1    0.00179  -0.00740 0.0116     9900. 0.704 Non-…
## 10 Clona… Plan… Simp… ClimPC.2   -0.00281  -0.0135  0.00890   10882. 0.629 Non-…
## 11 Clona… Plan… Simp… LHPC.1:LH…  0.000515 -0.00595 0.00681    9900  0.868 Non-…
## 12 Clona… Plan… Simp… ClimPC.1:… -0.00193  -0.00923 0.00569    9900. 0.619 Non-…
## 
## [[2]]
## # A tibble: 12 × 10
##    Trait     Taxa  Model Statistics post.mean    low95   high95 eff.samp   pMCMC
##    <chr>     <chr> <chr> <chr>          <dbl>    <dbl>    <dbl>    <dbl>   <dbl>
##  1 Growth_S… Plan… Phylo LHPC.1       2.18e-5 -2.99e-3  2.83e-3   10460. 9.80e-1
##  2 Growth_S… Plan… Phylo LHPC.2      -1.14e-3 -6.29e-3  3.79e-3    9900. 6.60e-1
##  3 Growth_S… Plan… Phylo ClimPC.1    -8.23e-4 -5.56e-3  4.03e-3    9900. 7.41e-1
##  4 Growth_S… Plan… Phylo ClimPC.2     5.96e-3  7.10e-4  1.13e-2    9900. 2.79e-2
##  5 Growth_S… Plan… Phylo LHPC.1:LH…  -9.03e-6 -2.55e-3  2.60e-3   10512. 9.85e-1
##  6 Growth_S… Plan… Phylo ClimPC.1:…  -3.48e-3 -6.81e-3 -1.62e-4    9900. 4.28e-2
##  7 Growth_S… Plan… Simp… LHPC.1      -2.10e-4 -1.62e-3  1.10e-3    9900. 7.66e-1
##  8 Growth_S… Plan… Simp… LHPC.2      -1.64e-3 -3.97e-3  1.01e-3    9900. 1.91e-1
##  9 Growth_S… Plan… Simp… ClimPC.1    -4.49e-4 -2.39e-3  1.39e-3    9900. 6.36e-1
## 10 Growth_S… Plan… Simp… ClimPC.2     4.79e-3  2.60e-3  7.09e-3    9900. 1.01e-4
## 11 Growth_S… Plan… Simp… LHPC.1:LH…  -2.11e-4 -1.39e-3  1.11e-3    9900. 7.41e-1
## 12 Growth_S… Plan… Simp… ClimPC.1:…  -2.78e-3 -4.27e-3 -1.31e-3    9900. 4.04e-4
## # ℹ 1 more variable: sig <chr>
## 
## [[3]]
## # A tibble: 12 × 10
##    Trait      Taxa  Model Statistics post.mean    low95  high95 eff.samp   pMCMC
##    <chr>      <chr> <chr> <chr>          <dbl>    <dbl>   <dbl>    <dbl>   <dbl>
##  1 Reproduct… Plan… Phylo LHPC.1     -0.00348  -1.09e-2 0.00413    9900. 0.350  
##  2 Reproduct… Plan… Phylo LHPC.2     -0.00360  -1.62e-2 0.00875    9900. 0.580  
##  3 Reproduct… Plan… Phylo ClimPC.1    0.00857  -1.20e-2 0.0297     9900. 0.415  
##  4 Reproduct… Plan… Phylo ClimPC.2    0.00166  -1.59e-2 0.0190    10617. 0.857  
##  5 Reproduct… Plan… Phylo LHPC.1:LH…  0.00425  -2.51e-3 0.0112     9900. 0.217  
##  6 Reproduct… Plan… Phylo ClimPC.1:…  0.0107   -1.69e-3 0.0226     9900. 0.0828 
##  7 Reproduct… Plan… Simp… LHPC.1     -0.000644 -1.21e-2 0.0111    10031. 0.912  
##  8 Reproduct… Plan… Simp… LHPC.2     -0.00214  -2.32e-2 0.0188     9594. 0.837  
##  9 Reproduct… Plan… Simp… ClimPC.1    0.0192    3.33e-3 0.0352     9900. 0.0222 
## 10 Reproduct… Plan… Simp… ClimPC.2   -0.00236  -2.07e-2 0.0173     9607. 0.801  
## 11 Reproduct… Plan… Simp… LHPC.1:LH…  0.0103    5.24e-6 0.0214     9186. 0.0610 
## 12 Reproduct… Plan… Simp… ClimPC.1:…  0.0168    3.86e-3 0.0290     9900. 0.00990
## # ℹ 1 more variable: sig <chr>
## 
## [[4]]
## # A tibble: 12 × 10
##    Trait  Taxa  Model Statistics post.mean    low95  high95 eff.samp pMCMC sig  
##    <chr>  <chr> <chr> <chr>          <dbl>    <dbl>   <dbl>    <dbl> <dbl> <chr>
##  1 Shrin… Plan… Phylo LHPC.1       7.20e-5 -2.63e-3 0.00297    9900. 0.958 Non-…
##  2 Shrin… Plan… Phylo LHPC.2      -7.43e-5 -5.04e-3 0.00463    9900. 0.977 Non-…
##  3 Shrin… Plan… Phylo ClimPC.1     7.47e-4 -4.09e-3 0.00528    8491. 0.738 Non-…
##  4 Shrin… Plan… Phylo ClimPC.2    -1.87e-4 -5.34e-3 0.00477    9900  0.949 Non-…
##  5 Shrin… Plan… Phylo LHPC.1:LH…   2.92e-4 -2.18e-3 0.00289    9900. 0.819 Non-…
##  6 Shrin… Plan… Phylo ClimPC.1:…  -6.99e-5 -3.39e-3 0.00316    9580. 0.966 Non-…
##  7 Shrin… Plan… Simp… LHPC.1       2.35e-4 -9.85e-4 0.00149    9900. 0.707 Non-…
##  8 Shrin… Plan… Simp… LHPC.2       1.28e-4 -2.15e-3 0.00239    9900. 0.913 Non-…
##  9 Shrin… Plan… Simp… ClimPC.1     5.21e-4 -1.25e-3 0.00221    9900. 0.546 Non-…
## 10 Shrin… Plan… Simp… ClimPC.2    -3.34e-4 -2.31e-3 0.00178    9900. 0.745 Non-…
## 11 Shrin… Plan… Simp… LHPC.1:LH…   3.15e-4 -8.41e-4 0.00150    9900. 0.595 Non-…
## 12 Shrin… Plan… Simp… ClimPC.1:…   1.12e-4 -1.28e-3 0.00145    9900. 0.871 Non-…
## 
## [[5]]
## # A tibble: 12 × 10
##    Trait     Taxa  Model Statistics post.mean    low95   high95 eff.samp   pMCMC
##    <chr>     <chr> <chr> <chr>          <dbl>    <dbl>    <dbl>    <dbl>   <dbl>
##  1 Survival… Plan… Phylo LHPC.1       -0.0145 -2.92e-2 -1.25e-4    9900  0.0511 
##  2 Survival… Plan… Phylo LHPC.2       -0.0176 -4.30e-2  8.64e-3    9900. 0.180  
##  3 Survival… Plan… Phylo ClimPC.1      0.0246  2.76e-3  4.62e-2    9900. 0.0277 
##  4 Survival… Plan… Phylo ClimPC.2      0.0181 -6.87e-3  4.17e-2    9900. 0.141  
##  5 Survival… Plan… Phylo LHPC.1:LH…    0.0135  9.09e-4  2.74e-2    9900. 0.0465 
##  6 Survival… Plan… Phylo ClimPC.1:…    0.0269  1.14e-2  4.33e-2    9900. 0.00101
##  7 Survival… Plan… Simp… LHPC.1       -0.0132 -2.65e-2  6.27e-4    9347. 0.0606 
##  8 Survival… Plan… Simp… LHPC.2       -0.0178 -4.23e-2  6.98e-3    9900. 0.165  
##  9 Survival… Plan… Simp… ClimPC.1      0.0249  5.83e-3  4.37e-2    9900  0.0101 
## 10 Survival… Plan… Simp… ClimPC.2      0.0171 -5.00e-3  3.87e-2    9596. 0.131  
## 11 Survival… Plan… Simp… LHPC.1:LH…    0.0136  1.99e-4  2.59e-2    9900. 0.0370 
## 12 Survival… Plan… Simp… ClimPC.1:…    0.0284  1.34e-2  4.27e-2    9900. 0.00101
## # ℹ 1 more variable: sig <chr>

11.2.2 Posterior distribution

Better for figures.

## Warning: Expected 1 pieces. Additional pieces discarded in 138600 rows [1, 2, 3, 4, 5,
## 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].

## Warning: Expected 1 pieces. Additional pieces discarded in 84 rows [1, 2, 3, 4, 5, 6, 7,
## 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].

11.3. Plotting the posterior distribution as ggplot2

Code

Code only

Cumulative

Figure 1

Vital rates

Figure 3

12. Estimation of Phylogenetic signal

12.1. Create an new ancillary function: `my.fake.lamb`

This function was created based on Garamszegi (2014)

12.2. Run the new function to estimate the the phylogenetic signal for all vital rates

The code below uses the my.fake.lamb to estimate phylogenetic signal for all models in a list.

The ggridges package is necessary (dependency for geom_density_ridges)

12.3. Comparison of the phylogenetic signal

without significant correlation between phylogenetic signal and contribution to buffering capacity

## 
##  Spearman's rank correlation rho
## 
## data:  filter(phylo_buffer_df, Taxa == "Plants")$MEDIAN and filter(phylo_buffer_df, Taxa == "Plants")$Relativemean
## S = 34, p-value = 0.2333
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##  rho 
## -0.7

Detailed phylogenetic signal for plants

Phylo x buffer plot

13. Testing the effect of timeserie length on buffering quantification

Effect size

## 
## Call:
## lm(formula = log(abs(Cumulative_SigElas)) ~ log(MatRep), data = final_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6043 -1.6915 -0.1736  1.7462  4.6299 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -5.1896     0.5375  -9.656   <2e-16 ***
## log(MatRep)   0.4773     0.3063   1.559    0.122    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.187 on 119 degrees of freedom
## Multiple R-squared:  0.02001,    Adjusted R-squared:  0.01177 
## F-statistic: 2.429 on 1 and 119 DF,  p-value: 0.1217

14. Final Metadata

Produces the final metadata table.

References

Bernard, Connor, Gabriel Silva Santos, Jacques A. Deere, Roberto Rodriguez-Caro, Pol Capdevila, Erik Kusch, Samuel J. L. Gascoigne, John Jackson, and Roberto Salguero-Gómez. 2023. “MOSAIC - A Unified Trait Database to Complement Structured Population Models.” Scientific Data 10 (1): 335. https://doi.org/10.1038/s41597-023-02070-w.

Demetrius, Lloyd. 1974. “Demographic Parameters and Natural Selection.” Proceedings of the National Academy of Sciences 71 (12): 4645–47. https://doi.org/10.1073/pnas.71.12.4645.

Garamszegi, László Zsolt. 2014. Modern Phylogenetic Comparative Methods and Their Application in Evolutionary Biology.

Jones, Owen R., Patrick Barks, Iain Stott, Tamora D. James, Sam Levin, William K. Petry, Pol Capdevila, et al. 2022. “Rcompadre and RageTwo R Packages to Facilitate the Use of the COMPADRE and COMADRE Databases and Calculation of Life-History Traits from Matrix Population Models.” Methods in Ecology and Evolution 13 (4): 770–81. https://doi.org/10.1111/2041-210X.13792.

Karger, Dirk Nikolaus, Olaf Conrad, Jürgen Böhner, Tobias Kawohl, Holger Kreft, Rodrigo Wilber Soria-Auza, Niklaus E. Zimmermann, H. Peter Linder, and Michael Kessler. 2017. “Climatologies at High Resolution for the Earth’s Land Surface Areas.” Scientific Data 4 (1): 170122. https://doi.org/10.1038/sdata.2017.122.

Note temperature is times 10. So, if temperature is 16ºC, data will show temperature value as 160↩︎
Remember that absolute sum of stochastic elasticity within respect to variance \(| \sum E_v^\sigma |\) is extremaly important to the analyses because our question is about the demographic position, not what determine the contribution of these vital rates. If we use raw values we could answers other question but not their position about buffering capacity↩︎

A general framework to run analyses in Santos et al.xxxxx

Gabriel Santos

2025-05-31 (20:10:36)

Readme

1 - Setting

0.1. SessionInfo and package versions

2. Data selection and cleaning - load population models & phylogenetic information

3. Life history traits (Calculate life history traits)

3.1 Fast-slow continuum / Life history PCA

4. Climate variables and environmental PCA

4.1. Extract and store climate information

4.1. Analyse climatic timeseries, calculate environmental stochasticity, and summarise information into a PCA

4.1.1. Checking collinearity and retaining non-redundant climatic variables

4.1.2. Climatic/environmental PCA

4.3 Mapping the axes of the environmental PCA

5. Calculate Demographic buffering from Vital Rates (Lower Level Parameters)

5.1. Sourcing main function used in this calculation

5.2. Derivative demographic buffering

5.2.1. Number of populations with timeseries longer than three

5.2.2. Merging Demographic Buffering into a single data frame

5.2.3. Produce a definitive demographic buffering data frame

6. Merging all (Climatic data + Life history + Buffering)

7. Building the super tree and making sure that works on phylogenetic analyses

8. Produce a final dataset and phylogenetic analyses

8.1. Prepare the associated phylogenetic information

9. GLMM parameters & data input

9.1. Prepare a separate dataset data_model

9.2. Selecting the best model

9.3. Determine the fixed effect of all models and response variables

9.4. OPTIMISING GLMM MODELS

10. Buffering patterns

10.1 Ranking buffering

Plot

10.3. Proportional contribution

10.3.1 Proportional contribution by vital rate

10.3.1 Proportional contribution by life form

10.3.3. Plots

Plot Elasticities

Plot proportional contribution

Raw figure

11. MCMCglmm results and tidy

11.1 Loading GLMM outputs

11.2. Data Harmonisation

11.2.1 As a summary

11.2.2 Posterior distribution

11.3. Plotting the posterior distribution as ggplot2

Code

Cumulative

Vital rates

12. Estimation of Phylogenetic signal

12.1. Create an new ancillary function: my.fake.lamb

12.2. Run the new function to estimate the the phylogenetic signal for all vital rates

12.3. Comparison of the phylogenetic signal

Detailed phylogenetic signal for plants

Phylo x buffer plot

13. Testing the effect of timeserie length on buffering quantification

Effect size

14. Final Metadata

References

5. Calculate Demographic buffering from Vital Rates
(Lower Level Parameters)

9.1. Prepare a separate dataset `data_model`

12.1. Create an new ancillary function: `my.fake.lamb`